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The  common  “folklore”  that  Eulerian-Lagrangian  methods  perform  better  (are  more  ac¬ 
curate)  with  large  Courant  numbers  (large  time  steps)  than  with  small  Courant  numbers, 
due  to  numerical  dispersion  in  the  latter  case,  is  explained  theoretically.  A  formulation 
that  does  not  suffer  from  large  numerical  dispersion  for  any  Courant  number  is  outlined. 

1.  INTRODUCTION 

There  appears  to  be  a  certain  amount  of  misunderstanding  in  the  hydrological  modeling 
community  of  the  extent  to  which  Eulerian-Lagrangian  methods  (ELMs)  for  advection- 
dispersion  transport  equations  are  numerically  dispersive.  One  often  hears  that  ELMs 
perform  well  for  problems  in  which  they  can  successfully  use  a  large  Courant  number 
(large  time  step),  but  that  when  forced  to  use  a  smaller  Courant  number,  and  hence  more 
time  steps,  they  suffer  from  numerical  dispersion  introduced  by  interpolation  at  each  step. 
Indeed,  some  formulations  of  ELMs  behave  in  precisely  this  way,  and  the  implementation 
and  use  of  such  formulations  has  probably  contributed  to  this  perception.  What  is  less 
well  understood  is  that  other  formulations  can  avoid  or  substantially  mitigate  this  effect. 

The  key  to  the  situation  is  the  relationship  between  averaging  in  the  mass  matrix  and 
averaging  in  the  right-hand-side  vector.  The  mass  matrix  represents  the  discrete  storage 
terms  at  the  new  time  level.  In  the  right-hand-side  vector,  advected  old-time  values  are 
averaged  by  interpolation  or  some  analogous  process  such  as  numerical  integration.  If 
the  averaging  on  these  two  sides  is  balanced,  then  fronts  can  be  propagated  with  minimal 
numerical  dispersion,  regardless  of  the  time  step.  In  numerically  dispersive  formulations, 
the  mass  matrix  is  typically  “lumped,”  i.e. ,  diagonal,  meaning  that  there  is  no  averaging 
in  the  mass  matrix  to  balance  the  averaging  on  the  right-hand  side. 

This  paper  will  detail  the  theoretical  basis  for  this  discussion,  outline  a  non-dispersive 
formulation,  and  discuss  the  expected  behavior  of  various  formulations  of  ELMs  in  some 
simple  cases. 

2.  A  SIMPLE  EULERIAN-LAGRANGIAN  METHOD 

The  concepts  to  be  discussed  in  this  paper  are  most  easily  and  clearly  presented  in  the 
context  of  a  1-D  constant-coefficient  advection-dispersion  equation,  where  there  are  no 
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extraneous  features.  Consider  the  model  problem 

Cu  =  ut  +  Vux  —  Duxx  =  uT  —  Duxx  =  f(x,  t ),  0  <  x  <  L,  t  >  0,  (1) 

with  appropriate  initial  and  boundary  conditions  that  are  not  relevant  here.  For  the  sake 
of  simplicity,  define  a  uniform  grid  with  Ax  =  L/I,  nodes  X*  =  iAx  ( i  =  0,1,...,/), 
and  cells  Cj  =  ^*+1/2]  =  [(*  —  1/2)  Ax,  ( i  +  1/2)  Ax]  centered  around  the  nodes. 

Similarly,  define  a  time  increment  At,  uniform  time  levels  tn  =  nAt  (n  =  0, 1,  2, . . .),  and 
time  intervals  Jn  =  [tn,tn+1]. 

The  space-time  r-directional  derivative  in  (1)  is  a  total  (or  Lagrangian,  substantial,  or 
material)  derivative  that  follows  the  advective  part  of  the  flow.  Along  this  direction,  with 
the  point  (x,tn+1)  we  associate  the  point  (x*,tn)  at  the  previous  time  level  that  flows  to 
(x,tn+1)  under  advection  (see  Figure  1),  namely 

x*  =  x  —  V  At.  (2) 


We  can  then  make  a  backward- Euler  approximation  of  the  total  derivative, 


uT (x, tn+1) 


u(x,  tn+1 )  —  u(x*,  tn ) 


Incorporating  (2)  and  (3)  into  a  centered  implicit  finite-difference  method  (FDM)  for  (1), 
in  the  interior  we  have  the  finite-difference  Eulerian-Lagrangian  method  (FDELM) 

ur+l-un(*t) ,  n-m1 + 2t/f+i  -  m1 


where  U  is  the  FDELM  approximation  to  u ,  f/”+1  =  U{xi,tn+l ),  and  f™+1  =  /(xj,tn+1). 
Let  5/t/j  =  Ui-i  —  2 Ui  +  Ui+i  denote  the  second  spatial  difference  operator  and  I  the 
identity  operator,  and  rewrite  (4)  in  the  form 


I  —  D- 


$)  ur1 = t'-w) + At/r1. 


Complete  specification  of  the  interior  FDELM  (5)  requires  evaluation  of  the  solution  at 
the  old  time  level,  Un(x*)  =  Un(xi  —  V At).  Because  x*  does  not  lie  at  a  node  in  general 
(see  Figure  1),  some  form  of  interpolation  or  the  equivalent  is  needed. 

Let  Cr  denote  the  Courant  number, 

Cr  =  VAt/Ax ,  (6) 

which  is  the  number  of  cells  traversed  by  advection  in  one  time  step.  For  definiteness, 
assume  that  Cr  >  0,  and  denote  its  integer  and  fractional  parts  by  [Cr]  and  (Cr), 
respectively.  The  simplest  evaluation  of  Un(x*)  is  by  linear  interpolation, 

Un(x*)  =  Un (xj  -  CrAx)  =  (1  -  (Cr))U?_[Cr]  +  (Cr)U?_[Cryi.  (7) 

Consider  the  case  0  <  Cr  <  1  ([Cr]  =  0,  (Cr)  =  Cr),  substitute  (7)  into  (5),  and  use  (6) 
to  obtain 


At  \  Un  —  TJn 

i  -  d^-sA  ur1  +  Atvu%  “'-1 

Ax2  XJ  *  Ax 


f/”  +  A  t/”+1; 
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in  (8)  we  have  precisely  a  first-order  explicit  upwind  FDM  in  the  advection  term.  This 
shows  that:  For  0  <  Cr  <  1,  the  FDELM  (5)  reduces  to  the  first-order  explicit  upwind 
FDM  (8).  This  method  is  known  to  be  numerically  dispersive. 

To  quantify  the  numerical  dispersivity  of  (8),  rewrite  it  in  the  form  of  a  second-order 
centered  FDM  plus  a  numerical  dispersion  term  [1]: 


D 


Tjn 

ur+i  +  Atv  *+l 


—  TIn 

ui-i 


2  Ax 


-  D  —  52Un 

^ num  FrV 


Ufi  +  a  tfl 


n+ 1 


(9) 


where  we  have  determined  the  numerical  dispersivity 


Dnum  =  v  Ax/ 2. 


(10) 


3.  DISCUSSION 

It  is  possible  for  the  FDELM  (5)  to  be  less  dispersive  and  more  accurate  with  larger  time 
steps,  Cr  >  1.  In  such  a  case,  in  a  time  step  the  solution  is  advected  [ Cr ]  cells  without 
dispersion,  then  is  evolved  by  the  upwind  FDM  to  complete  the  time  step.  Because  fewer 
interpolations  (upwind  steps)  are  needed  for  a  given  simulation  time  as  the  time  step  is 
increased,  the  resulting  U  is  less  dispersed.  This  effect  is  balanced  by  increasing  errors 
in  the  physical  dispersion  approximation  and  the  source  term  as  At  increases,  so  that  it 
is  plausible  to  seek  an  “optimal  time  step”  that  minimizes  the  sum  of  the  advective  and 
dispersive  errors  [2],  This  example,  and  others  similar  to  it,  appear  to  be  the  reason  for 
the  “folklore”  mentioned  in  the  abstract  and  introduction. 

Within  the  framework  of  (5),  numerical  dispersion  can  be  avoided  by  a  form  of  quadratic 
interpolation,  as  proved  in  [3].  There  a  second-order  error  estimate  was  demonstrated,  a 
consequence  of  which  is  that  first-order  numerical  dispersion  must  be  absent.  Rather  than 
pursuing  that  line  of  reasoning  further,  here  we  will  consider  finite- volume  ELMs  that  can 
alleviate  the  numerical  dispersion  in  a  systematic  way  that  is  relatively  straightforward 
to  understand. 

Neglecting  the  physical  dispersion  in  (5),  the  identity  operator  is  applied  to  Un+1,  and 
an  interpolation  operator  to  Un.  The  net  effect  is  that  the  solution  Un+1  will  be  some  type 
of  weighted  average  of  Un.  With  higher-order  interpolation  as  in  [3],  negative  weights  will 
be  necessary  if  numerical  dispersion  is  to  be  avoided,  creating  the  potential  for  spurious 
oscillations.  Nonlinear  interpolators  that  combine  linear  and  quadratic  techniques  or 
other  combinations,  such  as  slope  limiters,  can  be  used  [4],  These  are  more  difficult 
to  formulate  in  multiple  dimensions  and  near  boundaries.  In  ELMs  that  often  reach 
to  or  across  boundaries  in  determining  the  backtracked  point  x*  in  (2),  we  have  found 
it  natural  and  advantageous  to  modify  the  treatment  of  the  left-hand  side  of  (5),  rather 
than  the  right-hand  side.  This  will  be  described  in  the  context  of  a  finite-volume  Eulerian- 
Lagrangian  localized  adjoint  method  (FVELLAM)  [5,6]. 

4.  FINITE- VOLUME  ELLAM 

The  Eulerian-Lagrangian  localized  adjoint  method  (ELLAM)  [7]  uses  space-time  finite 
elements  oriented  along  the  r  directional  derivative  of  (1),  generalizing  ELMs  in  a  mass- 
conservative  way  that  systematically  treats  general  boundary  conditions.  This  is  one 


4 


X  X  X  X 

i-  512  i-312  i-112  i+112 


n+1 

t 


n 


t 


Figure  1.  An  FVELLAM  space-time  cell. 


of  many  possible  settings  for  the  ensuing  discussion,  but  seems  to  provide  some  insight 
on  possible  treatments  of  numerical  dispersion.  The  FVELLAM  uses  space-time  test 
functions  that  are  characteristic  (indicator)  functions  of  space-time  subdomains  oriented 
along  r.  This  will  be  described  briefly  here  to  set  the  stage. 

As  in  Figure  1,  define  the  cell  C*  =  [xj_i/2,  ^+1/2]  at  the  new  time  level  tn+1,  and  its 
backtracked  image  C*  =  [x*_1/2^x*+1/.2\  at  t,n,  and  the  space-time  cell  fi™+1  that  consists 
of  all  points  between  them.  A  discrete  interior  FVELLAM  equation  (we  do  not  consider 
boundary  issues  here)  is  obtained  by  integrating  (1)  over  fl”+1  in  space  and  time,  or 
equivalently  by  multiplying  (1)  by  the  characteristic  function  of  fl”+1  and  integrating 
over  the  whole  space-time  domain.  In  the  integration,  note  that  the  measures  dx  dt  and 
dx  dr  are  the  same,  so  that  we  can  think  of  “time”  integration  in  the  r  direction  instead 
of  t.  So  doing,  we  obtain 


uT  dr  dx  — 


Duxx  dx  dr  = 


(11) 


and  integrating  in  r  and  x,  respectively,  in  the  first  two  terms  of  (11), 


rXi- |_i/2  rtn+1 

/  un+1(x )  dx  +  D  [ux(xi- 1/2  -  V (tn+1  -  r),  r)  -  ux(xi+1/2  -  V (tn+1  -  r),  r)]  dr 

JXi- 1/2  Jt n 


f  1+1/2  un(x*)  dx*  +  [  f  dx  dr. 
Jx*  „„  V  ’  Jw+1 


(12) 


i- 1/2 


The  FVELLAM  is  obtained  by  choosing  approximations  Un  and  Un+l  for  un  and  un+1, 
and  numerical  integration  rules  in  (12).  We  use  a  backward-Euler  approximation  of  the 
time  integrals,  so  that  (12)  is  approximated  by 


PXi- {-1/2 

/  (7"+1(x)  dx  +  DAiKtyj+%  -  ([4)”+%] 

JXi- 1/2 

=  r<H/2  un(x*)  dx* + At  r+i/2  r+\x)  dx. 

J X*.  ,  •’  x  1  /o 


(13) 


The  crux  of  our  discussion  is  the  evaluation  of  the  Un+1  integral  on  the  left-hand  side  of 
(13)  and  the  Un  integral  on  the  right-hand  side. 
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Figure  2.  Area  showing  value  of  left-hand-side  Un+1  integral. 


4.1.  Un+1  integral 

In  a  finite-volume  method,  it  is  customary  to  think  of  the  solution  Un+1  as  piecewise- 
constant,  i.e.,  constant  on  each  spatial  cell  C*.  This  choice  would  lead  to  (5),  and  hence 
to  the  same  numerical  dispersion  noted  there.  Instead  we  view  Un+1  as  piecewise-linear, 
interpolated  between  nodes  x*.  The  Un+l  integral  in  (13)  is  then  the  area  depicted  in 
Figure  2,  which  is 

a; u"+iw ix = (\up' + \u"+i + \u^) = {r + isi) u"+i-  <i4> 

In  comparison  with  (5),  the  weights  1/8,  3/4,  1/8  in  (14)  “unlump”  the  left-hand  side 
of  (13).  If  we  approximate  Ux  by  a  difference  quotient  and  fn+1  by  a  piecewise-constant 
function,  we  can  write  (13)  in  the  form  analogous  to  (9), 

Ax  ('  -  (b  -  V"(x>)  dx’  +  AxA tfT1-  (15) 

We  see  that  unlumping  has  introduced  in  (15)  a  negative  numerical  dispersivity 
1  Ax2 

Dnum  ~  -g^p 

4.2.  Un  integral 

The  key  to  alleviating  numerical  dispersion  in  FVELLAM  is  to  balance  the  negative 
dispersivity  of  unlumping  with  positive  dispersivity  in  the  evaluation  of  the  Un  integral 
in  (13).  The  unlumping  has  created  “room”  for  a  weighted  average,  or  equivalently  a 
numerical  integration,  in  the  Un  integral  to  be  non-dispersive  in  the  end.  The  following 
discussion  is  directed  toward  achieving  the  desired  balance. 

For  the  simple  problem  considered  here,  exact  integration  is  possible,  evaluating  Un(x*) 
as  in  (7)  in  terms  of  the  Courant  number  Cr.  For  0  <  Cr  <  1/2,  C*  =  ^+1/2]  c 

[xi-i,  Xj+i],  so  that  a  single  formula  in  terms  of  U-l_1,  C7”,  and  U™+1  is  applicable.  A 
complementary  formula  covers  the  case  1/2  <  Cr  <  1,  and  simple  shifts  incorporate  the 
integer  part  [ Cr ].  For  0  <  Cr  <  1/2,  straightforward  algebra  yields 


L 


'<+ 1/2 


'1  Cr  Cr2' 


i- 1/2 


Un(x*)dx*  =  Axll-  +  — + 


Ui-i  +  [  j  ~  Cr2 


U? 
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+ 


Cr  Cr2\ 

“T  +  —  ) 


(17) 


For  Cr  =  0,  the  weights  1/8,  3/4,  1/8  reappear,  properly  balancing  the  left-hand  side.  In 
practice,  numerical  integration  will  be  necessary,  so  we  seek  a  quadrature  formula  that 
balances  well  over  a  range  of  Courant  numbers. 

To  conserve  global  mass  in  the  Un  integral,  in  general  it  is  necessary  to  evaluate  the 
integral  with  a  forward-tracking  procedure  [8,5].  Quadrature  points  are  regularly  dis¬ 
tributed  at  tn,  hence  integrate  the  total  mass  exactly.  Each  point  is  tracked  forward  to 
tn+1,  and  if  it  lands  in  C)  =  ^j+1/2],  its  mass  is  then  contributed  to  the  integral 

over  C*  =  [x\_1j2i  £*+1//2].  With  discontinuous  characteristic  functions  as  in  Figure  1,  the 
numerical  integral  can  be  changed  greatly  by  small  changes  in  Cr.  For  example,  with 
quadrature  points  Xi-i,  ^*—1/2,  2+  ^4-1/2,  £j+i  each  carrying  mass  (domain  length)  Ax/2, 
the  Un  integral  in  (13)  will  be  approximated  by  Ax(l/4  U/_1  +  3/4  U-1)  if  Cr  is  small 
and  positive,  but  by  Ax(3/4  UJ/  +  1/4  U-l+1)  if  Cr  is  small  and  negative.  This  discontin¬ 
uous  sensitivity  will  yield  unacceptable  results  in  practice,  so  for  numerical  purposes  it  is 
necessary  to  modify  the  integral  with  a  continuous  smoothed  test  function: 


f <+1/2  Un{x*)  dx*  «  [L  Un(x*)Wi(x* ,tn)  dx*, 
Jx*  ,  /0  JO 


i- 1/2 


(18) 


where  Wi(x*,tn)  is  a  continuous  approximation  of  the  discontinuous  characteristic  function 
of  C* .  A  natural  hrst  choice  is  to  take  Wi(x,tn+1)  to  be  the  piecewise-linear  hat  function 
equal  to  1  at  Xi  and  0  at  all  other  nodes,  and  let  Wi(x*,tn )  =  Wi(x,tn+1). 

The  piecewise-linear  14/  makes  the  integrand  on  the  right-hand  side  of  (18)  piecewise- 
quadratic.  If  integrated  exactly,  it  will  produce  weights  of  1/6,  2/3,  1/6  (the  weights 
of  Simpson’s  rule)  for  Cr  =  0.  These  do  not  match  the  desired  1/8,  3/4,  1/8.  In 
a  manner  analogous  to  (16),  a  right-hand-side  Dnum  =  1/6  Ax2 / At  results,  hence  an 
overall  numerical  dispersivity 


D 


num 


1\  Ax2 
8  )  ~ At ~ 


1  Ax2 
24  "At~' 


(19) 


For  fixed  Ax  and  small  At,  results  will  be  excessively  dispersed. 

Thus,  if  44/  is  used,  we  must  integrate  approximately  in  a  fashion  that  restores  the 
weights  1/8,  3/4,  1/8  of  the  exact  non- 14/  integral  for  the  case  Cr  =  0.  Comparing 
Figure  3  to  the  exact  non- 14/  area  in  Figure  2,  the  shaded  areas  are  the  same,  and  the 
shaded  area  in  Figure  3  would  be  obtained  from  the  trapezoidal  rule,  with  I4/(x*_i)  =  0, 
Wifa- 1/2)  =  1/2,  Wi(xi)  =  1,  Wi(xi+l/ 2)  =  1/2,  and  Wi(xi+1)  =  0.  For  0  <  Cr  <  1/2 
with  the  trapezoidal  rule,  we  have 


(L  Un(x*)WAx*,tn)dx* 
Jo 


-L 


A/((+H^-'+(3-f)c,"+G 

5+1/2 17-00  ar  +  (A  -  Af )  €u?, 


—^1  un 
4  J  Ul+1 


by  (17).  Analogous  to  (19),  the  overall  numerical  dispersivity  is 

1  1  2\  Ax2  V Ax  V2At 

-Cr - Cr 2  — —  = - . 

4  2  J  At  4  2 


D. 


(20) 


(21) 
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Figure  3.  Trapezoidal  rule  to  determine  correct  area,  Cr  =  0,  NS  =  2. 


Figure  4.  Trapezoidal  rule  to  determine  correct  area,  Cr  =  0,  NS  =  4. 


Unlike  (19),  as  At  — >  0  the  dispersivity  in  (21)  has  a  finite  limit  V  Ax/\1  which  is 
half  of  the  numerical  dispersivity  of  upwinding  in  (10).  Note  that  the  trapezoidal  rule  is 
calculated  on  a  mesh  of  size  Ax/ 2,  compared  to  upwinding  on  size  Ax  in  (8). 

Numerical  results  were  reported  in  Problem  1  of  [5]  for  a  test  case  with  V  —  25,  D  =  2.5, 
Ax  =  2,  At  =  0.001  (cell  Peclet  number  Pe  =  V Ax/D  —  20,  Cr  —  1/80),  L  =  250,  final 
time  2.  For  these  values,  (21)  gives  Dnum  =  12.2,  hence  the  total  dispersivity  of  14.7  is  5.9 
times  the  physical  dispersivity.  Then  the  numerical  front  width  should  be  about  a/579  ~ 
2.4  times  the  analytical  front  width.  The  widths  between  95%  and  5%  concentration  in 
the  numerical  and  analytical  solutions  were  lAAx  and  6Ax,  respectively.  For  the  same 
data  except  that  At  =  0.01  {Cr  =  1/8),  hence  Dnum  =  9.375,  D  +  Dnum  =  4.75D, 
\/4.75  ~  2.2,  the  numerical  and  analytical  front  widths  were  10A:r  and  5Ax,  respectively. 
This  indicates  that  the  analysis  explains  the  dispersivity  behavior  very  well. 

As  suggested  by  the  comparison  between  (21)  and  (10),  it  is  possible  to  reduce  the 
dispersivity  by  using  the  trapezoidal  rule  on  a  finer  set  of  quadrature  points.  Let  NS  be 
the  number  of  subintervals  in  the  cell  Cp  in  Figure  3,  NS  =  2.  Figure  4  shows  the  case 
NS  =  4.  The  exact  area  will  be  obtained  as  shown,  if  the  values  Wi{xi- 1)  =  IF* (£4-3/4)  = 
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Dr 


=  (  -Cr  —  -Cr2 
,8  2 


0,  Wi(xi- 1/2)  =  1/2,  =  Wi(xi)  =  Wi(xi+ 1/4)  =  1,  =  1/2,  and 

Wi(xi+ 3/4)  =  Wj(xi+i)  =  0  are  used  (this  is  depicted  in  Figure  4  of  [5]).  Thus,  instead  of 
a  hat  function,  W,t  in  this  case  changes  linearly  from  0  to  1  on  [xj_3/4,  x± -1/4],  and  from  1 
to  0  on  [xj+i/4,  Xj+3/4].  A  development  similar  to  (20)  concludes  that 

Ax2  UArr  V2At  UAz  ,  , 

— —  = - > - as  At  — >  0.  (22 

At  8  2  8  v  ' 

More  generally,  for  any  even  value  of  NS,  we  have  TT7* (^«— 1/2)  =  Wj(xj+i/2)  =  1/2, 

Wi  =  1  at  all  quadrature  points  in  between,  Wt  =  0  at  all  quadrature  points  outside,  and 
Wi  is  interpolated  linearly  between  quadrature  points.  With  these  specifications,  Dnum 
tends  to  the  limit  V Ax/2NS  as  At  — »  0.  Thus,  with  denser  trapezoidal-quadrature 
points,  numerical  dispersivity  can  be  made  very  small,  even  with  small  Courant  number. 


5.  CONCLUSIONS 

The  observed  numerical  dispersion  in  Eulerian-Lagrangian  methods  (ELMs),  particu¬ 
larly  with  small  time  steps,  is  not  an  intrinsic  feature  of  ELMs.  Rather,  it  is  a  result 
of  formulations  whose  mass  matrix  is  lumped,  as  is  typical  for  finite  differences.  With 
an  unlumped  mass  matrix,  and  careful  balancing  of  it  with  the  integration  rule  for  the 
old-time-level  mass,  one  can  design  ELMs  that  exhibit  minimal  numerical  dispersion  for 
any  Courant  number,  large  or  small. 
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